---
title: "Figure_Code_PP2025"
author: "Chloe Mortenson"
date: "2025-10-20"
output: html_document
---

################################################################################
# DEMOCRACY CONJOINT EXPERIMENT - VISUALIZATION SCRIPT
# 
# Purpose: Visualizations for all analyses including
#          Thurstone results, AMCE plots, interaction graphs, VIMP heatmaps,
#          IMCE density plots, and individual-level relationship graphs.
#
# Author: Chloe Mortenson & Dr.Erik Nisbet
# Date: 10/20/25
# Version: 1.0
#
# Prerequisites: Must run data analysis script first to generate model objects
#
# Input: Model objects and analysis results from data/ and models/ directories
# Output: High-resolution publication-ready figures in figures/ directory
#
################################################################################


```{r}
#===============================================================================
# ARTICLE FIGURES GENERATION SCRIPT
# This script generates all figures using outputs from the analysis script
# 
# PREREQUISITE: Run the complete analysis script first!
# That script creates all the required .rds files used here
#===============================================================================

# Load required libraries
library(dplyr)
library(ggplot2)
library(cregg)
library(cjbart)
library(tidyr)
library(gridExtra)
library(grid)
library(cowplot)
library(coda)
library(knitr)
library(kableExtra)

# Create figures directory if it doesn't exist
if (!dir.exists("figures")) {
  dir.create("figures")
}

cat("=== ARTICLE FIGURES GENERATION SCRIPT ===\n\n")

#===============================================================================
# LOAD ALL REQUIRED FILES
#===============================================================================

cat("Loading analysis outputs...\n\n")

# Load the data files created by the analysis script
model_data_fixed <- readRDS("final_fixed_data.rds")
cat("✅ Loaded: final_fixed_data.rds\n")

fixed_amce <- readRDS("final_fixed_amce.rds")
cat("✅ Loaded: final_fixed_amce.rds\n")

cjbart_model <- readRDS("final_fixed_cjbart.rds")
cat("✅ Loaded: final_fixed_cjbart.rds\n")

fixed_imce <- readRDS("final_fixed_imce.rds")
cat("✅ Loaded: final_fixed_imce.rds\n")

vimp_results <- readRDS("final_vimp_results.rds")
cat("✅ Loaded: final_vimp_results.rds\n")

cat("\n✅ All data loaded successfully!\n\n")

# Display summary information
cat("Data summary:\n")
cat("  - Observations:", nrow(model_data_fixed), "\n")
cat("  - Respondents:", length(unique(model_data_fixed$ResponseId)), "\n")
cat("  - IMCE respondents:", nrow(fixed_imce$imce), "\n")

#===============================================================================
# VERIFY ALL REQUIRED FILES EXIST
#===============================================================================

cat("Checking for required data files from analysis script...\n\n")

required_files <- c(
  "final_fixed_data.rds",
  "final_fixed_amce.rds",
  "final_fixed_cjbart.rds",
  "final_fixed_imce.rds",
  "final_vimp_results.rds"
)

missing_files <- c()
for(file in required_files) {
  if(file.exists(file)) {
    cat("✅", file, "\n")
  } else {
    cat("❌", file, "MISSING\n")
    missing_files <- c(missing_files, file)
  }
}

if(length(missing_files) > 0) {
  cat("\n⚠️ ERROR: Missing required files!\n")
  cat("Please run the complete analysis script first to generate:\n")
  for(file in missing_files) {
    cat("  -", file, "\n")
  }
  stop("Cannot proceed without required data files.")
}

cat("\n✅ All required files present.\n\n")

#===============================================================================
# PREPARE VIMP DATA FOR VISUALIZATION
#===============================================================================

cat("Processing VIMP data for visualization...\n")

# Debug: Check the structure of vimp_results
cat("Checking VIMP structure...\n")
str(vimp_results)

# Extract VIMP results from the nested structure
vimp_df <- as.data.frame(vimp_results$results$results)

cat("Successfully extracted VIMP data frame with columns:", paste(names(vimp_df), collapse=", "), "\n")

# Verify we have the expected columns
expected_cols <- c("Attribute", "Level", "covar", "importance")
missing_cols <- setdiff(expected_cols, names(vimp_df))
if(length(missing_cols) > 0) {
  stop(paste("Missing required columns:", paste(missing_cols, collapse=", ")))
}

# Process VIMP results for visualization
working_vimp_data <- vimp_df %>%
  mutate(
    importance = pmin(importance * 100, 100),
    lab_text = sprintf("%.1f", importance),
    covar_label = case_when(
      covar == "Liberal_Score_const" ~ "Liberal Democracy",
      covar == "Procedural_Score_const" ~ "Procedural Democracy", 
      covar == "Distributive_Score_const" ~ "Distributive Democracy",
      covar == "educ_const" ~ "Education",
      covar == "polint_const" ~ "Political Interest",
      covar == "ideo_const" ~ "Ideology",
      covar == "age_const" ~ "Age",
      TRUE ~ covar
    ),
    Level_Clean = case_when(
      grepl("Majority", Level) ~ "Majority",
      grepl("Minority", Level) ~ "Minority",
      grepl("Normative", Level) ~ "Normative",
      TRUE ~ Level
    )
  ) %>%
  filter(Level_Clean %in% c("Majority", "Minority")) %>%
  filter(
    !is.na(importance) & 
    !is.infinite(importance) & 
    importance >= 0 & 
    !is.na(Attribute) & 
    !is.na(covar_label)
  ) %>%
  mutate(
    importance = pmin(importance, 100),
    lab_text = sprintf("%.1f", importance),
    Attribute_Display = case_when(
      Attribute == "Presidential_Role" ~ "Political Equality",
      Attribute == "Press_Freedom" ~ "Freedom of Expression", 
      Attribute == "Rule_Law" ~ "Rule of Law",
      Attribute == "Econ_Wellbeing" ~ "Economic Well-being",
      TRUE ~ Attribute
    )
  )

# Create factor levels for proper ordering
working_vimp_data$covar_label <- factor(
  working_vimp_data$covar_label, 
  levels = c(
    "Liberal Democracy", "Procedural Democracy", "Distributive Democracy",
    "Ideology", "Age", "Education", "Political Interest"
  )
)

cat("✅ VIMP data prepared\n")
cat("  - VIMP variables:", nrow(working_vimp_data), "\n\n")

#===============================================================================
# SECTION 1: DIAGNOSTIC CHECKS
#===============================================================================

cat("SECTION 1: DIAGNOSTIC CHECKS\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("1.1 Creating CJBART convergence diagnostics...\n")

# Check sigma convergence
if (!is.null(cjbart_model$sigma)) {
  sigma_mcmc <- mcmc(cjbart_model$sigma)
  geweke_result <- geweke.diag(sigma_mcmc)
  
  cat("Geweke Z-score for sigma:", geweke_result$z, "\n")
  if (abs(geweke_result$z) < 1.96) {
    cat("✅ Chain has converged (|Z| < 1.96)\n")
  } else {
    cat("⚠️  Chain may not have converged (|Z| > 1.96)\n")
  }
  
  # Create diagnostic plots
  png("figures/diagnostic_sigma_trace.png", width = 10, height = 8, units = "in", res = 300)
  par(mfrow=c(2,1))
  plot(sigma_mcmc, main="Trace plot of sigma parameter")
  acf(cjbart_model$sigma, main="Autocorrelation of sigma parameter")
  par(mfrow=c(1,1))
  dev.off()
  
  cat("✅ Saved: figures/diagnostic_sigma_trace.png\n")
}

# Check varcount convergence
if (!is.null(cjbart_model$varcount)) {
  varcount_sums <- rowSums(cjbart_model$varcount)
  varcount_mcmc <- mcmc(varcount_sums)
  geweke_result <- geweke.diag(varcount_mcmc)
  
  cat("Geweke Z-score for variable usage:", geweke_result$z, "\n")
  
  png("figures/diagnostic_varcount_trace.png", width = 10, height = 8, units = "in", res = 300)
  par(mfrow=c(2,1))
  plot(varcount_mcmc, main="Trace plot of variable usage")
  acf(varcount_sums, main="Autocorrelation of variable usage")
  par(mfrow=c(1,1))
  dev.off()
  
  cat("✅ Saved: figures/diagnostic_varcount_trace.png\n")
}

cat("\n")

#===============================================================================
# PREPARE DATA FOR MARGINAL MEANS PLOTS
#===============================================================================

cat("Preparing data for marginal means plots...\n")

# Prepare data for MM plots
model_data_clean <- model_data_fixed %>%
  mutate(
    Econ_Wellbeing = factor(Econ_Wellbeing, levels = c("Econ_Minority", "Econ_Normative", "Econ_Majority")),
    Press_Freedom = factor(Press_Freedom, levels = c("Press_Minority", "Press_Normative", "Press_Majority")),
    Rule_Law = factor(Rule_Law, levels = c("Law_Minority", "Law_Normative", "Law_Majority")),
    Presidential_Role = factor(Presidential_Role, levels = c("Pres_Minority", "Pres_Normative", "Pres_Majority")),
    selected = as.numeric(selected)
  )

cat("✅ Data prepared for marginal means plots\n\n")

# Function to create MM plots
create_mm_plot <- function(mm_data, title, legend_title, by_var_label = "Economic Wellbeing") {
  # Clean both BY and level columns to remove prefixes
  mm_data_clean <- mm_data %>%
    mutate(
      BY_clean = case_when(
        grepl("Minority", BY) ~ "Minority",
        grepl("Normative", BY) ~ "Normative",
        grepl("Majority", BY) ~ "Majority",
        TRUE ~ BY
      ),
      BY_clean = factor(BY_clean, levels = c("Minority", "Normative", "Majority")),
      level_clean = case_when(
        grepl("Minority", level) ~ "Minority",
        grepl("Normative", level) ~ "Normative",
        grepl("Majority", level) ~ "Majority",
        TRUE ~ level
      ),
      level_clean = factor(level_clean, levels = c("Majority", "Minority", "Normative")),
      # Create position adjustment - Minority centered (0), Normative left (-0.25), Majority right (+0.25)
      position_offset = case_when(
        level_clean == "Minority" ~ 0,
        level_clean == "Normative" ~ -0.25,
        level_clean == "Majority" ~ 0.25
      )
    )
  
  ggplot(mm_data_clean, aes(x = BY_clean, y = estimate, color = level_clean, group = level_clean)) +
    geom_point(aes(x = as.numeric(BY_clean) + position_offset), size = 2.5) +
    geom_errorbar(aes(x = as.numeric(BY_clean) + position_offset, ymin = lower, ymax = upper),
                 width = 0.2) +
    theme_minimal() +
    scale_color_manual(
      values = c("Minority" = "#FC8D62", "Normative" = "#8DA0CB", "Majority" = "#66C2A5"),
      breaks = c("Majority", "Minority", "Normative"),
      name = legend_title
    ) +
    scale_x_continuous(
      breaks = 1:3,
      labels = c("Minority", "Normative", "Majority")
    ) +
    labs(title = title, y = "Probability of Profile Selection", x = by_var_label) +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5, size = 12),
      axis.text = element_text(size = 9),
      axis.title = element_text(size = 10),
      legend.position = "bottom"
    )
}
  

# Function to safely save and display plots
safe_plot_save <- function(plot_obj, filename, width = 8, height = 5) {
  tryCatch({
    ggsave(filename, plot_obj, width = width, height = height, dpi = 300)
    cat("✅ Saved:", filename, "\n")
    print(plot_obj)
  }, error = function(e) {
    cat("⚠️  Error saving", filename, ":", conditionMessage(e), "\n")
    cat("Attempting alternative save method...\n")
    tryCatch({
      png(filename, width = width * 300, height = height * 300, res = 300)
      print(plot_obj)
      dev.off()
      cat("✅ Saved using alternative method:", filename, "\n")
    }, error = function(e2) {
      cat("❌ Could not save plot:", conditionMessage(e2), "\n")
    })
  })
}

#===============================================================================
# FIGURE 2: MARGINAL MEANS OF ALL DEMOCRATIC ATTRIBUTES
#===============================================================================

cat("FIGURE 2: MARGINAL MEANS OF ALL DEMOCRATIC ATTRIBUTES\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 2: Marginal means for all attributes...\n")

# Calculate marginal means for all attributes
mm_econ <- cj(data = model_data_clean, formula = selected ~ Econ_Wellbeing,
             id = ~ResponseId, estimate = "mm")
mm_press <- cj(data = model_data_clean, formula = selected ~ Press_Freedom,
              id = ~ResponseId, estimate = "mm")
mm_law <- cj(data = model_data_clean, formula = selected ~ Rule_Law,
            id = ~ResponseId, estimate = "mm")
mm_pres <- cj(data = model_data_clean, formula = selected ~ Presidential_Role,
             id = ~ResponseId, estimate = "mm")

# Add attribute labels
mm_econ$attribute <- "Economic Wellbeing"
mm_press$attribute <- "Freedom of Expression"
mm_law$attribute <- "Rule of Law"
mm_pres$attribute <- "Political Equality"

# Combine all marginal means
all_mm <- bind_rows(mm_econ, mm_press, mm_law, mm_pres) %>%
  mutate(
    level_clean = case_when(
      grepl("Minority", level) ~ "Minority",
      grepl("Normative", level) ~ "Normative",
      grepl("Majority", level) ~ "Majority",
      TRUE ~ level
    ),
    level_clean = factor(level_clean, levels = c("Minority", "Normative", "Majority")),
    attribute = factor(attribute, levels = c("Economic Wellbeing", "Freedom of Expression", 
                                              "Political Equality", "Rule of Law"))
  )

# Create combined marginal means plot matching the style in the image
fig2_all_mm <- ggplot(all_mm, aes(x = level_clean, y = estimate, color = attribute, group = attribute)) +
  geom_point(size = 3, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
               width = 0.15, position = position_dodge(width = 0.4), linewidth = 0.8) +
  theme_minimal() +
  scale_color_manual(
    values = c("Economic Wellbeing" = "#E41A1C", 
               "Freedom of Expression" = "#377EB8",
               "Political Equality" = "#4DAF4A",
               "Rule of Law" = "#984EA3"),
    name = "Attribute"
  ) +
  labs(
    title = "Marginal Means for Attributes of Democracy",
    x = "Level",
    y = "Probability of Choice"
  ) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 14),
    axis.text = element_text(size = 11, color = "black"),
    axis.title = element_text(size = 12, face = "bold"),
    legend.position = "right",
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA)
  ) +
  scale_y_continuous(limits = c(0.40, 0.61), breaks = seq(0.40, 0.60, by = 0.05))

safe_plot_save(fig2_all_mm, "figures/fig2_all_attributes_mm.png", width = 10, height = 7)

cat("\n")

#===============================================================================
# FIGURE 3: RULE OF LAW BY ECONOMIC WELLBEING
#===============================================================================

cat("FIGURE 3: RULE OF LAW BY ECONOMIC WELLBEING\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 3: Rule of Law by Economic Wellbeing...\n")

law_mm <- cj(data = model_data_clean, formula = selected ~ Rule_Law,
            id = ~ResponseId, estimate = "mm", by = ~Econ_Wellbeing)
law_plot <- create_mm_plot(law_mm, "Rule of Law by Economic Wellbeing", "Levels of Rule of Law")
safe_plot_save(law_plot, "figures/fig3_rule_of_law_interaction.png")

cat("\n")

#===============================================================================
# FIGURE 4: POLITICAL EQUALITY BY ECONOMIC WELLBEING
#===============================================================================

cat("FIGURE 4: POLITICAL EQUALITY BY ECONOMIC WELLBEING\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 4: Political Equality by Economic Wellbeing...\n")

pres_mm <- cj(data = model_data_clean, formula = selected ~ Presidential_Role,
             id = ~ResponseId, estimate = "mm", by = ~Econ_Wellbeing)
pres_plot <- create_mm_plot(pres_mm, "Political Equality by Economic Wellbeing", "Levels of Political Equality")
safe_plot_save(pres_plot, "figures/fig4_political_equality_interaction.png")

cat("\n")

#===============================================================================
# FIGURE 5: FREEDOM OF EXPRESSION BY ECONOMIC WELLBEING
#===============================================================================

cat("FIGURE 5: FREEDOM OF EXPRESSION BY ECONOMIC WELLBEING\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 5: Freedom of Expression by Economic Wellbeing...\n")

press_mm <- cj(data = model_data_clean, formula = selected ~ Press_Freedom,
              id = ~ResponseId, estimate = "mm", by = ~Econ_Wellbeing)
press_plot <- create_mm_plot(press_mm, "Freedom of Expression by Economic Wellbeing", "Levels of Freedom of Expression")
safe_plot_save(press_plot, "figures/fig5_press_freedom_interaction.png")

cat("\n")

#===============================================================================
# FIGURE 6: IMCE DENSITY PLOTS
#===============================================================================

cat("FIGURE 6: IMCE DENSITY PLOTS (REVEALED PREFERENCES)\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 6: IMCE density plots...\n")

# Reshape IMCE data
imces_data <- fixed_imce$imce %>%
  select(ResponseId, Press_Minority, Press_Majority, Pres_Majority, Pres_Minority,
         Law_Minority, Law_Majority, Econ_Majority, Econ_Minority) %>%
  pivot_longer(cols = -ResponseId, names_to = "attribute", values_to = "imce") %>%
  mutate(
    attribute_group = case_when(
      grepl("Press_", attribute) ~ "Freedom of Expression",
      grepl("Law_", attribute) ~ "Rule of Law",
      grepl("Econ_", attribute) ~ "Economic Well-being",
      grepl("Pres_", attribute) ~ "Political Equality"
    ),
    level = ifelse(grepl("Majority", attribute), "Majority", "Minority")
  )

# Function to create density plots
make_density_plot <- function(data, title) {
  ggplot(data, aes(x = imce, fill = level, color = level)) +
    geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
    geom_density(alpha = 0.3) +
    theme_minimal() +
    labs(title = title, x = "", y = "Density") +
    scale_x_continuous(limits = c(-0.4, 0.1)) +
    theme(plot.title = element_text(size = 12, hjust = 0.5),
          legend.position = "right",
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 8)) +
    scale_fill_manual(values = c("red", "blue"), name = "Level") +
    scale_color_manual(values = c("red", "blue"), name = "Level")
}

# Create individual plots
p1 <- make_density_plot(imces_data %>% filter(attribute_group == "Freedom of Expression"), 
                       "Freedom of Expression")
p2 <- make_density_plot(imces_data %>% filter(attribute_group == "Rule of Law"), "Rule of Law")
p3 <- make_density_plot(imces_data %>% filter(attribute_group == "Economic Well-being"), 
                       "Economic Well-being")
p4 <- make_density_plot(imces_data %>% filter(attribute_group == "Political Equality"), 
                       "Political Equality")

# Combine and save
combined_densities <- arrangeGrob(p1, p2, p3, p4, ncol = 1)
ggsave("figures/fig6_imce_densities.png", combined_densities, width = 8, height = 10, dpi = 300)

cat("✅ Saved: figures/fig6_imce_densities.png\n\n")

#===============================================================================
# FIGURE 7: VIMP HEATMAP
#===============================================================================

cat("FIGURE 7: IMPORTANCE SCORE INTERACTION HEATMAP\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 7: VIMP heatmap...\n")

vimp_heatmap <- ggplot(working_vimp_data, 
                       aes(x = covar_label, y = Level_Clean, fill = importance)) +
  facet_grid(Attribute_Display ~ ., space = "free", scales = "free", switch = "y") +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = lab_text), 
            color = ifelse(working_vimp_data$importance > 50, "white", "black"),
            size = 4, fontface = "bold") +
  scale_fill_gradient(low = "white", high = "firebrick2", 
                     breaks = c(25, 50, 75, 100), limits = c(0, 100),
                     name = "Importance (%)") +
  scale_y_discrete(limits = c("Minority", "Majority")) +
  labs(x = "Democracy Understanding and Individual Attributes",
       y = "Attribute Level",
       title = "Variable Importance in Predicting Individual Conjoint Preferences") +
  theme_minimal() +
  theme(
    text = element_text(size = 14),
    plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14, face = "bold"),
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
    panel.grid = element_blank(),
    strip.text.y.left = element_text(angle = 0, hjust = 1, size = 12, face = "bold"),
    strip.background = element_rect(fill = "grey90", color = NA),
    panel.spacing = unit(1, "lines"),
    strip.placement = "outside"
  )

ggsave("figures/fig7_vimp_heatmap.png", vimp_heatmap, width = 12, height = 8, dpi = 300)
print(vimp_heatmap)
cat("✅ Saved: figures/fig7_vimp_heatmap.png\n\n")

#===============================================================================
# INDIVIDUAL-LEVEL RELATIONSHIP PLOTS FUNCTION
#===============================================================================

# Function to create quartile relationship graphs
create_quartile_graph <- function(imce_data, outcome_var, covariate_var, 
                                 covariate_label, outcome_label, output_filename) {
  
  plot_data <- tibble(
    ResponseId = imce_data$imce$ResponseId,
    outcome = imce_data$imce[[outcome_var]],
    covariate = imce_data$imce[[covariate_var]]
  ) %>% filter(!is.na(covariate) & !is.na(outcome))
  
  quartiles <- quantile(plot_data$covariate, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
  
  plot_data <- plot_data %>%
    mutate(
      covar_quartile = case_when(
        covariate <= quartiles[2] ~ "Q1 (Lowest)",
        covariate <= quartiles[3] ~ "Q2", 
        covariate <= quartiles[4] ~ "Q3",
        TRUE ~ "Q4 (Highest)"
      ),
      covar_quartile = factor(covar_quartile, 
                             levels = c("Q1 (Lowest)", "Q2", "Q3", "Q4 (Highest)"))
    ) %>%
    arrange(outcome) %>%
    mutate(effect_order = 1:n())
  
  quartile_counts <- plot_data %>%
    count(covar_quartile) %>%
    mutate(percentage = n / sum(n) * 100,
           label = paste0("(n=", n, ", ", sprintf("%.1f", percentage), "%)"))
  
  quartile_labels <- setNames(quartile_counts$label, quartile_counts$covar_quartile)
  
  overall_se <- sd(plot_data$outcome, na.rm = TRUE) / sqrt(nrow(plot_data))
  
  plot_data <- plot_data %>%
    mutate(smooth_mean = mean(outcome),
           lower = smooth_mean - 1.96 * overall_se,
           upper = smooth_mean + 1.96 * overall_se)
  
  y_min <- floor(min(plot_data$outcome, na.rm = TRUE) * 10) / 10
  y_max <- ceiling(max(plot_data$outcome, na.rm = TRUE) * 10) / 10
  
  # Top panel: Effect line
  effect_line <- ggplot(plot_data, aes(x = effect_order, y = outcome)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey70", alpha = 0.5) +
    geom_line(color = "black", linewidth = 0.8) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    labs(x = "", y = outcome_label,
         title = paste("Individual Marginal Component Effects:", gsub("Effect", "", outcome_label)),
         subtitle = paste("Reference Group: Normative", gsub(" \\(.*\\) Effect", "", outcome_label))) +
    theme_minimal() +
    theme(panel.grid = element_blank(),
          plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
          plot.subtitle = element_text(size = 12, hjust = 0.5),
          axis.text.y = element_text(size = 11),
          axis.title.y = element_text(angle = 90, vjust = 0.5, size = 12),
          axis.text.x = element_blank(), 
          axis.ticks.x = element_blank()) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(limits = c(y_min, y_max), 
                      breaks = pretty(c(y_min, y_max), n = 6),
                      labels = function(x) sprintf("%.2f", x))
  
  # Bottom panel: Quartile histogram
  quartile_hist <- ggplot(plot_data, aes(x = effect_order)) +
    geom_histogram(aes(fill = covar_quartile), bins = 100, position = "stack", boundary = 0) +
    labs(x = "Individual Respondents", y = "Count",
         fill = paste(covariate_label, "Quartiles")) +
    scale_fill_manual(
      values = c("Q1 (Lowest)" = "#66C2A5", "Q2" = "#FC8D62", 
                "Q3" = "#8DA0CB", "Q4 (Highest)" = "#E78AC3"),
      labels = function(x) paste0(x, "\n", quartile_labels[x])
    ) +
    theme_minimal() +
    theme(panel.grid = element_blank(), 
          legend.position = "bottom",
          legend.box.margin = margin(t = 10),
          axis.text.y = element_text(size = 11),
          axis.title = element_text(size = 12),
          legend.title = element_text(size = 11),
          legend.text = element_text(size = 9)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    guides(fill = guide_legend(nrow = 1))
  
  # Combine panels
  combined_plot <- plot_grid(effect_line, quartile_hist, ncol = 1, 
                            rel_heights = c(0.6, 0.4), align = "v")
  
  ggsave(output_filename, combined_plot, dpi = 300, width = 12, height = 8)
  
  return(combined_plot)
}

#===============================================================================
# FIGURE 8: LIBERAL DEMOCRACY SCORE & ECONOMIC WELLBEING (MAJORITY)
#===============================================================================

cat("FIGURE 8: LIBERAL DEMOCRACY & ECONOMIC WELLBEING (MAJORITY)\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 8: Liberal Democracy Score & Economic Wellbeing (Majority)...\n")

p8 <- create_quartile_graph(
  imce_data = fixed_imce,
  outcome_var = "Econ_Majority",
  covariate_var = "Liberal_Score_const",
  covariate_label = "Liberal Democracy Score", 
  outcome_label = "Economic Well-being (Majority) Effect",
  output_filename = "figures/fig8_liberal_econ_majority.png"
)

cat("✅ Saved: figures/fig8_liberal_econ_majority.png\n\n")

#===============================================================================
# FIGURE 9: LIBERAL DEMOCRACY SCORE & ECONOMIC WELLBEING (MINORITY)
#===============================================================================

cat("FIGURE 9: LIBERAL DEMOCRACY & ECONOMIC WELLBEING (MINORITY)\n")
cat("-" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("Creating Figure 9: Liberal Democracy Score & Economic Wellbeing (Minority)...\n")

p9 <- create_quartile_graph(
  imce_data = fixed_imce,
  outcome_var = "Econ_Minority",
  covariate_var = "Liberal_Score_const",
  covariate_label = "Liberal Democracy Score",
  outcome_label = "Economic Well-being (Minority) Effect", 
  output_filename = "figures/fig9_liberal_econ_minority.png"
)

cat("✅ Saved: figures/fig9_liberal_econ_minority.png\n\n")

#===============================================================================
# SUMMARY
#===============================================================================

cat("\n")
cat("=" %>% rep(80) %>% paste(collapse = ""), "\n")
cat("FIGURE GENERATION COMPLETE\n")
cat("=" %>% rep(80) %>% paste(collapse = ""), "\n\n")

cat("All figures saved to the 'figures/' directory:\n\n")

cat("DIAGNOSTIC FIGURES:\n")
cat("  - diagnostic_sigma_trace.png\n")
cat("  - diagnostic_varcount_trace.png\n\n")

cat("MAIN ANALYSIS FIGURES:\n")
cat("  Figure 2: fig2_all_attributes_mm.png (Marginal Means of All Democratic Attributes)\n")
cat("  Figure 3: fig3_rule_of_law_interaction.png (Rule of Law by Economic Wellbeing)\n")
cat("  Figure 4: fig4_political_equality_interaction.png (Political Equality by Economic Wellbeing)\n")
cat("  Figure 5: fig5_press_freedom_interaction.png (Freedom of Expression by Economic Wellbeing)\n")
cat("  Figure 6: fig6_imce_densities.png (IMCE Revealed Preferences for Democracy)\n")
cat("  Figure 7: fig7_vimp_heatmap.png (Importance Score Interaction Heatmap)\n")
cat("  Figure 8: fig8_liberal_econ_majority.png (Liberal Democracy & Econ Wellbeing - Majority)\n")
cat("  Figure 9: fig9_liberal_econ_minority.png (Liberal Democracy & Econ Wellbeing - Minority)\n\n")

cat("✅ All article figures ready for publication!\n")
```



